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Abstract. We present a new multi-fluid, grid MHD code PIERNIK, 
which is based on the Relaxing TVD scheme ( |Jin and Xin, 1995| . The 
original scheme (see Trac & Pen ()2003p and Pen et al. ((2003])) has been 
extended by an addition of dynamically independent, but interacting 
fluids: dust and a diffusive cosmic ray gas, described within the fluid 
approximation, with an option to add other fluids in an easy way. 
The code has been equipped with shearing-box boundary conditions, 
and a selfgravity module, Ohmic resistivity module, as well as other 
facilities which are useful in astrophysical fluid-dynamical simulations. 
The code is parallelized by means of the MPI library. In this paper 
we introduce the multifluid extension of Relaxing TVD scheme and 
present a test case of dust migration in a two-fluid disk composed of 
gas and dust. We demonstrate that due to the difference in azimuthal 
velocities of gas and dust and the drag force acting on both components 
dust drifts towards maxima of gas pressure distribution. 



1 Multifluid extension of the Relaxing TVD scheme 

The basic set of conservative MHD equations (see Paper I, ( [Hanasz et al., 2008^ , 
this volume) describes a single fluid. The Relaxing TVD scheme by Pen, Arras 
& Wong l|2003p can be easily extended for multiple fluids by concatenation of the 
vectors of conservative variables for different fluids 
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representing ionized gas, neutral gas, as well as dust treated as a pressureless fluid. 
In a short notation this can be written as 

u= (u^u",u''). (1.2) 

The corresponding fluxes for these fluids are combined in a similar way 

F(u,B) = (F''(u',B),F"(u"),F''(u'')) , (1.3) 

where the elementary flux vectors like F*(u*, B), F"(u") and F'^(u'') are considered 
as fluxes computed independently for each fluid in the single fluid description. In 
multidimensional computations the fluxes G(u, B) and H(u, B), corresponding to 
the transport of conservative quantities in y and z-dircctions, are constructed in 
a similar manner. The full multifluid system of MHD equations, including source 
terms 

atu + a,F = s(u), (1.4) 

is subsequently solved by means of the Relaxing TVD scheme described in Paper I, 
together with the induction equation coupling magnetic filed B to the ionized gas 
component. The term S(u) includes any source terms (like gravity) corresponding 
to individual fluids as well as the terms coupling the dynamics of different fluids. 

As an example of mutual fluid interaction we consider the environment of 
protoplanetary discs, where aerodynamic interaction of gas and dust particles takes 
place. The interaction of the neutral gas and dust components induces the effects 
of the drag force in the momentum and energy equations for these fluids 

St = d^-p'p- - ^) , = ^St, (1.5) 

5," = a-^^Vt^n • - ^) , (1-6) 

where 

(1-7) 

is the inverse product of gas density and particle stopping time. In general a''" 
depends on properties of dust grains (theirs sizes and bulk densities) and gas, and 
their relative velocities. For simplicity in the calculations presented we assumed its 
constant value of 10. Friction forces can be incorporated, when needed, between 
any pair of fluids. 



2 Dust migration in protoplanetary disks 

Radial migration of dust particles through gaseous protoplanetary disks is essential 
for planetesimals formation (see e.g. Johansen et al. (|2007p and references therein). 
In this section we present a test example of dust migration in a 2D gaseous disk, 
under the action of a drag force, coupling the dust and gas components, defined in 
formula (jl.Sp . We assume that the two-component disk rotates in the gravitational 
field of a point mass. Initially gas forms a bell-like distribution, with a density 
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Fig. 1. Left panel: Initial density radial distribution of dust (solid line) and gas (dashed 
line). Right panel: Initial azimuthal velocity of dust (solid line) and gas (dashed line). 
Dust velocity is Keplerian, outside the softening radius of point mass gravity. Gas velocity 
is super-keplerian inside radius 0.4 and sub-keplerian outside due to the gas pressure 
gradient. 
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Fig. 2. Initial gas density (left panel) and dust density (right panel) at t = 0. A quarter 
of full disk disk is simulated with the aid of the "corner-periodic" boundary conditions. 

maximum in the mid of the disk radius, while the dust component is uniformly 
distributed across the disk (see Figs. [1] and [5]). The simulation is performed in 
a Cartesian domain with a spatial resolution Ux x Uy — 256 x 256. To speed up 
the simulation we use 'corner-periodic' boundary conditions, which enforce 90° 
periodicity of simulated objects in the azimuthal angle. 

In order to reduce numerical artifacts near the outer domain corner we force 
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gas, outside the disk radius, to follow the assumed rotation curve through addi- 
tional drag terms added in each timestep to x- and y-components of momentum 

A(puj) = -adamp(wi -'yj,mit)pAt, i^{x,y} (2.1) 

where Vi and Wijnit are the current and the initial velocity components respectively, 
Q^damp S [0, 1] is a factor responsible for artificial damping velocity fluctuations 
(deviations from the initial velocity distribution) on disk's outskirts. 

A radial gradient of gas density in a protoplanetary disks leads to a sub- 
keplerian rotation. On the other hand dust treated as pressureless fluid tends 
to orbit the central mass with the Kcplcrian velocity. As a result, there is a 
difference between the azimuthal velocities of gas and dust in regions where the 
radial gradient of gas pressure is non- vanishing. Since dust interacts with gas by 
means of the friction force, one can expect an exchange of momenta between gas 
and dust, and thus dust migration in the radial direction ( [Rice et al., 2004| ). In 
the present disc configuration the gas rotation is faster than Keplerian inside the 
radius of gas density maximum and slower outside, due to the pressure gradient 
contribution to the radial force balance. In the inner region gas speeds up the 
dust rotation and in the outer region the dust rotation is slowed down, thus the 
resulting torques shift dust towards the gas pressure maximum. 




Fig. 3. Left panel: Dust density distribution after 1.5 rotation periods at the radius of 
the initial maximum gas pressure. The dust component apparently gathers at the radius 
corresponding to the gas pressure maximum (where density gradient is zero). Right 
panel: Line plots of radial distribution of dust (solid line) and gas (dashed line). Values 
for dust are magnified 5 times for plotting on the right panel. 

The results of our experiment, performed with the PIERNIK code, are shown in 
Fig. El The effect of dust migration towards the gas pressure maximum is apparent, 
as expected, since after two rotational periods the majority of dust accumulates in 
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a ring of large density. The simulation presented in this paper demonstrates that 
the Relaxing TVD scheme can be successfully extended to describe the dynamics 
of multiple fluids interacting through the drag force. 
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